! include 'mkl_pardiso.f90'

module MySparseOperator
    implicit none
    integer, parameter :: mso_dp = 8    ! double presion



!     interface mso_sparse_transpose
!         module procedure mso_sparse_d_transpose, mso_sparse_z_transpose
!     end interface ! mso_sparse_transpose


contains



    subroutine mso_sparse_d_create_coo(rowind, colind, acoo, A, dimA, cols)
        !*******************************************************************************
        !  Content:
        !  mso_sparse_d_create_coo: create a double sparse matrix from sparse matrix in COO form
        !*******************************************************************************
        use common_functions, only: cf_qsort_i, cf_allocate, cf_deallocate
        use ifport
        implicit none
        integer, dimension(:), intent(in) :: rowind, colind
        real(kind=8), dimension(:), intent(in) :: acoo
        type(mso_dsparse_matrix), intent(inout) :: A
        integer, intent(in) :: dimA  ! dimensions of A
        integer, intent(in), optional :: cols

        integer :: nnz, job(8), info, idump
        integer, dimension(:), allocatable :: indr, pst
        integer :: ii, jj
        real(kind=8), dimension(:), allocatable :: dtmp
        real(kind=8) ddump



        if ( size(rowind).ne.size(colind) .or. size(colind).ne.size(acoo) ) then
            write(*,*) 'the size of rowind, colind, acoo should be agree with each other in mso_sparse_d_create_coo!'
            stop 1
        end if
        if ( minval(rowind).lt.1 .or. maxval(rowind).gt.dimA ) then
            write(*,*) 'rowind in mso_sparse_d_create_coo is not defined well in mso_sparse_d_create_coo!'
            stop 1
        end if
        if ( .not.present(cols) ) then
            if ( minval(colind).lt.1 .or. maxval(colind).gt.dimA ) then
                write(*,*) 'colind in mso_sparse_d_create_coo is not defined well in mso_sparse_d_create_coo!'
                stop 1
            end if
        else
            if ( minval(colind).lt.1 .or. maxval(colind).gt.cols ) then
                write(*,*) 'colind in mso_sparse_d_create_coo is not defined well in mso_sparse_d_create_coo!'
                stop 1
            end if
        end if

        A%rows = dimA
        if ( .not.present(cols) ) then
            A%cols = dimA
        else
            A%cols = cols
        end if

        !         if ( .not.present(cols) ) then
        !             A%cols = dimA
        !         else
        !             A%cols = cols
        !         end if

        !         nnz = size(rowind)

        !         call cf_allocate(dtmp,nnz)
        !         dtmp = acoo
        !         call cf_allocate(indr,nnz)
        !         call cf_allocate(pst,nnz)
        !         indr = (/ 1:nnz /)
        !         pst = (rowind-1)*A%cols + colind

        !         !         call cf_qsort_i(indr,nnz,md_cmp)
        !         call qsort_i4(indr,nnz,kind(1),md_cmp)

        !         ii = 1
        !         jj = 2
        !         do while ( jj.le.nnz )
        !           if ( pst(indr(ii)).eq.pst(indr(jj)) ) then
        !               dtmp(indr(ii)) = dtmp(indr(ii)) + dtmp(indr(jj))
        !               jj = jj + 1
        !           else
        !               if ( jj.gt.ii+1 ) then
        !                   ii = ii + 1
        !                   indr(ii) = indr(jj)
        !                   jj = jj + 1
        !               else
        !                   ii = ii + 1
        !                   jj = jj + 1
        !               end if
        !           end if
        !         end do
        !         nnz = ii

        !         call cf_allocate(A%a,nnz)
        !         call cf_allocate(A%ja,nnz)
        !         call cf_allocate(A%ia,A%rows+1)

        !         A%a = dtmp(indr(1:nnz))
        !         A%ja = colind(indr(1:nnz))

        !         A%ia = 0
        !         do ii = 1, nnz, 1
        !             A%ia(rowind(indr(ii))+1) = A%ia(rowind(indr(ii))+1) + 1
        !         end do
        !         A%ia(1) = 1
        !         do ii = 2, A%rows+1, 1
        !             A%ia(ii) = A%ia(ii-1) + A%ia(ii)
        !         end do
        job(1) = 2  ! converted to CSR and he column indices in CSR representation are sorted in the increasing order within each row.
        job(2) = 1  ! one-based indexing for CSR
        job(3) = 1  ! one-based indexing for COO
        job(4) = 2
        job(6) = 3
        call cf_allocate(A%ia,A%rows+1)
        nnz = size(rowind)
        if ( .not.allocated(A%ja) .or. (size(A%ja).lt.nnz) ) then
            call cf_allocate(A%ja,nnz)
        end if
        if ( .not.allocated(A%a) .or. (size(A%a).lt.nnz) ) then
            call cf_allocate(A%a,nnz)
        end if
        call mkl_dcsrcoo( job, A%rows, A%a, A%ja, A%ia, nnz, acoo, rowind, colind, info)


        !         call cf_deallocate(indr)
        !         call cf_deallocate(pst)
        !         call cf_deallocate(dtmp)
        !     contains
        !       function md_cmp(a,b)
        !           integer, intent(in) :: a,b
        !           integer(2) :: md_cmp
        !           if ( pst(a).lt.pst(b) ) md_cmp = -1
        !           if ( pst(a).eq.pst(b) ) md_cmp = 0
        !           if ( pst(a).gt.pst(b) ) md_cmp = 1
        !       end function md_cmp
    end subroutine mso_sparse_d_create_coo

    subroutine mso_sparse_z_create_coo(rowind, colind, acoo, A, dimA, cols)
        !*******************************************************************************
        !  Content:
        !  mso_sparse_z_create_coo: create a double complex sparse matrix from sparse matrix in COO form
        !*******************************************************************************
        use common_functions, only: cf_qsort_i, cf_allocate, cf_deallocate
        implicit none
        integer, dimension(:), intent(in) :: rowind, colind
        complex(kind=8), dimension(:), intent(in) :: acoo
        type(mso_zsparse_matrix), intent(inout) :: A
        integer, intent(in) :: dimA  ! dimensions of A
        integer, intent(in), optional :: cols

        integer :: nnz, job(8), info, err
        integer, dimension(:), allocatable :: indr, pst
        integer :: ii, jj
        complex(kind=8), dimension(:), allocatable :: dtmp


        if ( size(rowind).ne.size(colind) .or. size(colind).ne.size(acoo) ) then
            write(*,*) 'the size of rowind, colind, acoo should be agree with each other in mso_sparse_z_create_coo!'
            stop 1
        end if
        if ( minval(rowind).lt.1 .or. maxval(rowind).gt.dimA ) then
            write(*,*) 'rowind in mso_sparse_d_create_coo is not defined well in mso_sparse_z_create_coo!'
            stop 1
        end if
        if ( .not.present(cols) ) then
            if ( minval(colind).lt.1 .or. maxval(colind).gt.dimA ) then
                write(*,*) 'colind in mso_sparse_z_create_coo is not defined well in mso_sparse_z_create_coo!'
                stop 1
            end if
        else
            if ( minval(colind).lt.1 .or. maxval(colind).gt.cols ) then
                write(*,*) 'colind in mso_sparse_z_create_coo is not defined well in mso_sparse_z_create_coo!'
                stop 1
            end if
        end if

        A%rows = dimA
        if ( .not.present(cols) ) then
            A%cols = dimA
        else
            A%cols = cols
        end if

        nnz = size(rowind)

        call cf_allocate(dtmp,nnz)
        dtmp = acoo
        call cf_allocate(indr,nnz)
        call cf_allocate(pst,nnz)
        indr = (/ 1:nnz /)
        pst = (rowind-1)*A%cols + colind

        call cf_qsort_i(indr,nnz,md_cmp)

        ii = 1
        jj = 2
        do while ( jj.le.nnz )
            if ( pst(indr(ii)).eq.pst(indr(jj)) ) then
                dtmp(indr(ii)) = dtmp(indr(ii)) + dtmp(indr(jj))
                jj = jj + 1
            else
                if ( jj.gt.ii+1 ) then
                    ii = ii + 1
                    indr(ii) = indr(jj)
                    jj = jj + 1
                else
                    ii = ii + 1
                    jj = jj + 1
                end if
            end if
        end do
        nnz = ii

        call cf_allocate(A%a,nnz)
        call cf_allocate(A%ja,nnz)
        call cf_allocate(A%ia,A%rows+1)

        A%a = dtmp(indr(1:nnz))
        A%ja = colind(indr(1:nnz))

        A%ia = 0
        do ii = 1, nnz, 1
            A%ia(rowind(indr(ii))+1) = A%ia(rowind(indr(ii))+1) + 1
        end do
        A%ia(1) = 1
        do ii = 2, A%rows+1, 1
            A%ia(ii) = A%ia(ii-1) + A%ia(ii)
        end do


        !         job(1) = 2  ! converted to CSR and he column indices in CSR representation are sorted in the increasing order within each row.
        !         job(2) = 1  ! one-based indexing for CSR
        !         job(3) = 1  ! one-based indexing for COO
        !         job(6) = 3

        !         call mkl_zcsrcoo( job, dimA, A%a, A%ja, A%ia, nnz, dtmp(indr(1:nnz)), rowind(indr(1:nnz)), colind(indr(1:nnz)), info)


        call cf_deallocate(indr)
        call cf_deallocate(pst)
        call cf_deallocate(dtmp)
    contains
        function md_cmp(a,b)
            integer, intent(in) :: a,b
            integer :: md_cmp
            if ( pst(a).lt.pst(b) ) md_cmp = -1
            if ( pst(a).eq.pst(b) ) md_cmp = 0
            if ( pst(a).gt.pst(b) ) md_cmp = 1
        end function md_cmp
    end subroutine mso_sparse_z_create_coo

    subroutine mso_sparse_d_transpose(A,B,trans)
        ! trans: 'T', 't', 'C', 'c'
        use common_functions, only: cf_allocate, cf_deallocate
        implicit none
        type(mso_dsparse_matrix), intent(in) :: A
        type(mso_dsparse_matrix), intent(inout) :: B
        character(len=1), intent(in), optional :: trans

        integer :: ii, jj, nnz
        integer, dimension(:), allocatable :: rowind

        if ( present(trans) .and. trans.ne.'t' .and. trans.ne.'T' .and. trans.ne.'C' .and. trans.ne.'c' ) then
            write(*,*) 'trans is not well defined in mso_sparse_d_transpose!'
            stop 1
        end if

        nnz = A%ia(A%rows+1) - 1
        call cf_allocate(rowind,nnz)
        do ii = 1, A%rows, 1
            do jj = A%ia(ii), A%ia(ii+1)-1, 1
                rowind(jj) = ii
            end do
        end do

        call mso_sparse_d_create_coo(A%ja(1:nnz),rowind,A%a(1:nnz),B,A%cols, A%rows)

        call cf_deallocate(rowind)
    end subroutine mso_sparse_d_transpose

    subroutine mso_sparse_z_transpose(A,B,trans)
        ! trans: 'T', 't', 'C', 'c'
        use common_functions, only: cf_allocate, cf_deallocate
        implicit none
        type(mso_zsparse_matrix), intent(in) :: A
        type(mso_zsparse_matrix), intent(inout) :: B
        character(len=1), intent(in), optional :: trans

        integer :: ii, jj
        integer, dimension(:), allocatable :: rowind

        if ( present(trans) .and. trans.ne.'t' .and. trans.ne.'T' .and. trans.ne.'C' .and. trans.ne.'c' ) then
            write(*,*) 'trans is not well defined in mso_sparse_z_transpose!'
            stop 1
        end if

        call cf_allocate(rowind,size(A%a))
        do ii = 1, A%rows, 1
            do jj = A%ia(ii), A%ia(ii+1)-1, 1
                rowind(jj) = ii
            end do
        end do

        if ( .not.present(trans) .or. trans.eq.'T' .or. trans.eq.'t' ) then
            call mso_sparse_z_create_coo(A%ja,rowind,A%a,B,A%cols,A%rows)
            return
        end if

        call mso_sparse_z_create_coo(A%ja,rowind,dconjg(A%a),B,A%cols,A%rows)

        call cf_deallocate(rowind)
    end subroutine mso_sparse_z_transpose


    subroutine mso_sparse_d_setype(A,mtype)
        ! 1 real and structurally symmetric
        ! 2 real and symmetric positive definite
        ! -2 real and symmetric indefinite
        ! 11 real and nonsymmetric
        type(mso_dsparse_matrix), intent(inout) :: A
        integer, intent(in) :: mtype

        A%mtype = mtype
    end subroutine mso_sparse_d_setype

    subroutine mso_sparse_z_setype(A,mtype)
        ! 3 complex and structurally symmetric
        ! 4 complex and Hermitian positive definite
        ! -4 complex and Hermitian indefinite
        ! 6 complex and symmetric
        ! 13 complex and nonsymmetric
        type(mso_zsparse_matrix), intent(inout) :: A
        integer, intent(in) :: mtype

        A%mtype = mtype
    end subroutine mso_sparse_z_setype

    subroutine mso_sparse_d_copy(A,B)
        use common_functions, only: cf_allocate
        type(mso_dsparse_matrix), intent(in) :: A
        type(mso_dsparse_matrix), intent(inout) :: B

        integer :: nnz

        if ( A%rows*A%cols.eq.0 ) then
            call mso_sparse_d_destroy(B)
            return
        end if

        nnz = A%ia(A%rows+1) - 1

        if ( .not.allocated(B%a) .or. (size(B%a).lt.nnz) ) then
            call cf_allocate(B%a,nnz)
        end if
        B%a(1:nnz) = A%a(1:nnz)

        if ( .not.allocated(B%ja) .or. (size(B%ja).lt.nnz) ) then
            call cf_allocate(B%ja,nnz)
        end if
        B%ja(1:nnz) = A%ja(1:nnz)

        if ( .not.allocated(B%ia) .or. (size(B%ia).ne.(A%rows+1)) ) then
            call cf_allocate(B%ia,A%rows+1)
        end if
        B%ia = A%ia

        B%rows = A%rows
        B%cols = A%cols
        B%mtype = A%mtype
    end subroutine mso_sparse_d_copy

    subroutine mso_sparse_z_copy(A,B)
        type(mso_zsparse_matrix), intent(in) :: A
        type(mso_zsparse_matrix), intent(inout) :: B

        integer :: nnz

        if ( A%rows*A%cols.eq.0 ) then
            call mso_sparse_z_destroy(B)
            return
        end if

        nnz = size(A%a)

        if ( allocated(B%a) ) then
            if ( size(B%a).ne.nnz ) then
                deallocate(B%a)
                allocate(B%a(nnz))
            end if
        else
            allocate(B%a(nnz))
        end if
        B%a = A%a

        if ( allocated(B%ja) ) then
            if ( size(B%ja).ne.nnz ) then
                deallocate(B%ja)
                allocate(B%ja(nnz))
            end if
        else
            allocate(B%ja(nnz))
        end if
        B%ja = A%ja

        if ( allocated(B%ia) ) then
            if ( size(B%ia).ne.size(A%ia) ) then
                deallocate(B%ia)
                allocate(B%ia(A%rows+1))
            end if
        else
            allocate(B%ia(A%rows+1))
        end if
        B%ia = A%ia

        B%rows = A%rows
        B%cols = A%cols
        B%mtype = A%mtype
    end subroutine mso_sparse_z_copy

    subroutine mso_sparse_d_dm(beta,A,B)
        ! B = beta*A, beta is double
        type(mso_dsparse_matrix), intent(in) :: A
        type(mso_dsparse_matrix), intent(inout) :: B
        real(kind=8), intent(in) :: beta

        integer :: nnz

        if ( A%rows*A%cols.eq.0 ) then
            call mso_sparse_d_destroy(B)
            return
        end if

        call mso_sparse_d_copy(A,B)
        B%a = beta*B%a
    end subroutine mso_sparse_d_dm

    subroutine mso_sparse_z_zm(beta,A,B)
        ! B = beta*A, beta is double complex
        type(mso_zsparse_matrix), intent(in) :: A
        type(mso_zsparse_matrix), intent(inout) :: B
        complex(kind=8), intent(in) :: beta

        if ( A%rows*A%cols.eq.0 ) then
            call mso_sparse_z_destroy(B)
            return
        end if

        call mso_sparse_z_copy(A,B)
        B%a = beta*B%a
    end subroutine mso_sparse_z_zm

    subroutine mso_sparse_d_add(A,beta,B,C,op)
        ! C = A + beta*op(B), op = 'N' or 'n' or 'T' or 't'
        ! op = 'N' or 'n' means C=A+beta*B
        ! op = 'T' or 't' means C=A+beta*B^T
        use common_functions, only: cf_allocate
        type(mso_dsparse_matrix), intent(in) :: A, B
        type(mso_dsparse_matrix), intent(inout) :: C
        real(kind=8), intent(in) :: beta
        character, intent(in), optional :: op

        character(len=1) :: trans
        integer :: nzmax, info, ii

        if ( present(op) ) then
            trans = op
        else
            trans = 'N'
        end if

        if ( A%rows.eq.0 .or. A%cols.eq.0 ) then
            if ( B%rows.eq.0 .or. B%cols.eq.0 ) then
                call mso_sparse_destroy(C)
            else
                call mso_sparse_copy(B,C)
                C%a = beta * C%a
            end if
            return
        end if

        if ( B%rows.eq.0 .or. B%cols.eq.0 ) then
            call mso_sparse_copy(A,C)
            return
        end if

        if ( trans.eq.'N' .or. trans.eq.'n' ) then
            if ( A%rows.ne.B%rows .or. A%cols.ne.B%cols ) then
                write(*,*) 'mso_sparse_d_add: the sizes of A and B should agree with each other!'
                stop 1
            end if
        else
            if ( A%rows.ne.B%cols .or. A%cols.ne.B%rows ) then
                write(*,*) 'mso_sparse_d_add: the sizes of A and B should agree with each other!'
                stop 1
            end if
        end if

        C%rows = A%rows
        C%cols = A%cols
        if ( .not.allocated(C%ia) .or. (size(C%ia).ne.A%rows+1) ) then
            call cf_allocate(C%ia,A%rows+1)
        end if

        call mkl_dcsradd(trans,1,0,A%rows,A%cols, A%a,A%ja,A%ia, beta, B%a,B%ja,B%ia, C%a,C%ja,C%ia,nzmax, info)

        nzmax = C%ia(A%rows+1) - 1
       if ( .not.allocated(C%a) .or. (size(C%a).lt.nzmax) ) then
            call cf_allocate(C%a, nzmax)
        end if
        if ( .not.allocated(C%ja) .or. (size(C%ja).lt.nzmax) ) then
            call cf_allocate(C%ja, nzmax)
        end if
        call mkl_dcsradd(trans,2,0,A%rows,A%cols, A%a,A%ja,A%ia, beta, B%a,B%ja,B%ia, C%a,C%ja,C%ia,nzmax, info)
    end subroutine mso_sparse_d_add

    subroutine mso_sparse_z_add(A,beta,B,C,op)
        ! C = A + beta*op(B), op = 'N' or 'n' or 'C' or 'c'
        ! op = 'N' or 'n' means C=A+beta*B
        ! op = 'T' or 't' means C=A+beta*B^T
        ! op = 'C' or 'c' means C=A+beta*B^H
        type(mso_zsparse_matrix), intent(in) :: A, B
        type(mso_zsparse_matrix), intent(inout) :: C
        complex(kind=8), intent(in) :: beta
        character, intent(in), optional :: op

        character(len=1) :: trans
        integer :: nzmax, info, ii

        if ( present(op) ) then
            trans = op
        else
            trans = 'N'
        end if

        if ( A%rows.eq.0 .or. A%cols.eq.0 ) then
            if ( B%rows.eq.0 .or. B%cols.eq.0 ) then
                call mso_sparse_destroy(C)
            else
                call mso_sparse_copy(B,C)
                C%a = beta * C%a
            end if
            return
        end if

        if ( B%rows.eq.0 .or. B%cols.eq.0 ) then
            call mso_sparse_copy(A,C)
            return
        end if

        if ( trans.eq.'N' .or. trans.eq.'n' ) then
            if ( A%rows.ne.B%rows .or. A%cols.ne.B%cols ) then
                write(*,*) 'mso_sparse_z_add: the sizes of A and B should agree with each other!'
                stop 1
            end if
        else
            if ( A%rows.ne.B%cols .or. A%cols.ne.B%rows ) then
                write(*,*) 'mso_sparse_z_add: the sizes of A and B should agree with each other!'
                stop 1
            end if
        end if

        C%rows = A%rows
        C%cols = A%cols

        if ( allocated(C%ia) ) then
            if ( size(C%ia).ne.A%rows+1 ) then
                deallocate(C%ia)
                allocate(C%ia(A%rows+1))
            end if
        else
            allocate(C%ia(A%rows+1))
        end if
        call mkl_zcsradd(trans,1,0,A%rows,A%cols, A%a,A%ja,A%ia, beta, B%a,B%ja,B%ia, C%a,C%ja,C%ia,nzmax, info)

        nzmax = C%ia(A%rows+1) - 1
        if ( allocated(C%a) ) then
            if ( size(C%a).ne.nzmax ) then
                deallocate(C%a)
                allocate(C%a(nzmax))
            end if
        else
            allocate(C%a(nzmax))
        end if
        if ( allocated(C%ja) ) then
            if ( size(C%ja).ne.nzmax ) then
                deallocate(C%ja)
                allocate(C%ja(nzmax))
            end if
        else
            allocate(C%ja(nzmax))
        end if
        call mkl_zcsradd(trans,2,0,A%rows,A%cols, A%a,A%ja,A%ia, beta, B%a,B%ja,B%ia, C%a,C%ja,C%ia,nzmax, info)
    end subroutine mso_sparse_z_add

    subroutine mso_sparse_d_add2(A,beta,B,op)
        ! A = A + beta*op(B), op = 'N' or 'n' or 'T' or 't'
        ! op = 'N' or 'n' means A=A+beta*B
        ! op = 'T' or 't' means A=A+beta*B^T
        implicit none
        type(mso_dsparse_matrix), intent(inout) :: A
        real(kind=8) :: beta
        type(mso_dsparse_matrix), intent(in) :: B
        character, intent(in), optional :: op

        type(mso_dsparse_matrix) :: C
        character :: trans

        if ( present(op) ) then
            trans = op
        else
            trans = 'N'
        end if
        call mso_sparse_d_add(A,beta,B,C,trans)
        call mso_sparse_d_copy(C,A)


        call mso_sparse_d_destroy(C)
    end subroutine mso_sparse_d_add2

    subroutine mso_sparse_z_add2(A,beta,B,op)
        ! A = A + beta*op(B), op = 'N' or 'n' or 'T' or 't'
        ! op = 'N' or 'n' means A=A+beta*B
        ! op = 'T' or 't' means A=A+beta*B^T
        implicit none
        type(mso_zsparse_matrix), intent(inout) :: A
        complex(kind=8) :: beta
        type(mso_zsparse_matrix), intent(in) :: B
        character, intent(in), optional :: op

        type(mso_zsparse_matrix) :: C
        character :: trans

        if ( present(op) ) then
            trans = op
        else
            trans = 'N'
        end if

        call mso_sparse_z_add(A,beta,B,C,trans)
        call mso_sparse_z_copy(C,A)

        call mso_sparse_destroy(C)
    end subroutine mso_sparse_z_add2

    subroutine mso_sparse_d_deleteRowCol(A,v,RC,B)
        type(mso_dsparse_matrix), intent(in) :: A
        integer, dimension(:), intent(in) :: v
        character, intent(in) :: RC
        type(mso_dsparse_matrix), intent(inout) :: B

        integer :: n, ii, sv, m, deletedcols

        call mso_dcsrdltRC(A%a,A%ia,A%ja,v,RC, B%a,B%ia,B%ja, deletedcols)
        B%rows = size(B%ia) - 1
        B%cols = A%cols - deletedcols
    end subroutine mso_sparse_d_deleteRowCol

    subroutine mso_sparse_z_deleteRowCol(A,v,RC,B)
        type(mso_zsparse_matrix), intent(in) :: A
        integer, dimension(:), intent(in) :: v
        character, intent(in) :: RC
        type(mso_zsparse_matrix), intent(inout) :: B

        integer :: n, ii, sv, m, deletedcols

        call mso_zcsrdltRC(A%a,A%ia,A%ja,v,RC, B%a,B%ia,B%ja, deletedcols)
        B%rows = size(B%ia) - 1
        B%cols = A%cols - deletedcols
    end subroutine mso_sparse_z_deleteRowCol

    subroutine mso_sparse_d_mv(A,x,y,trans)
        ! y = trans(A)*x
        use common_functions, only: cf_allocate, cf_deallocate
        type(mso_dsparse_matrix), intent(in) :: A
        real(kind=8), dimension(:), intent(in) :: x
        real(kind=8), dimension(:), allocatable, intent(inout) :: y
        character, optional :: trans

        character :: op
        integer, dimension(:), allocatable :: rowind
        integer :: nnz, ii, jj

        if ( .not.present(trans) ) then
            op = 'N'
        else
            op = trans
        end if

        if ( op.eq.'N' .or. op.eq.'n' ) then
            if ( size(x).ne.A%cols ) then
                write(*,*) 'size(x) should be equal to A%cols in mso_sparse_d_mv!'
                stop 1
            end if
        else
            if ( op.eq.'t' .or. op.eq.'T' .or. op.eq.'C' .or. op.eq.'c' ) then
                if ( size(x).ne.A%rows ) then
                    write(*,*) 'size(x) should be equal to A%rows in mso_sparse_d_mv when trans is set to be N or n!'
                    stop 1
                end if
            else
                write(*,*) 'trans is not well defined in mso_sparse_d_mv!'
                stop 1
            end if
        end if

        if ( op.eq.'N' .or. op.eq.'n' ) then
            call cf_allocate(y,A%rows)
        else
            call cf_allocate(y,A%cols)
        end if

        call mkl_dcsrgemv(op,A%rows,A%a,A%ia,A%ja,x,y)


        !         y = 0.d0

        !         if ( op.eq.'N' .or. op.eq.'n' ) then
        !             do ii = 1, A%rows, 1
        !                 do jj = A%ia(ii), A%ia(ii+1)-1, 1
        !                     y(ii) = y(ii) + A%a(jj)*x(A%ja(jj))
        !                 end do
        !             end do
        !             return
        !         end if

        !         nnz = size(A%a)

        !         call cf_allocate(rowind,nnz)
        !         do ii = 1, A%rows, 1
        !             do jj = A%ia(ii), A%ia(ii+1)-1, 1
        !                 rowind(jj) = ii
        !             end do
        !         end do

        !         do ii = 1, nnz, 1
        !             y(A%ja(ii)) = y(A%ja(ii)) + A%a(ii)*x(rowind(ii))
        !         end do
        !         call cf_deallocate(rowind)
    end subroutine mso_sparse_d_mv

    subroutine mso_sparse_z_mv(A,x,y,trans)
        ! y = trans(A)*x, A, x, y are all double complex
        use common_functions, only: cf_allocate, cf_deallocate
        type(mso_zsparse_matrix), intent(in) :: A
        complex(kind=8), dimension(:), intent(in) :: x
        complex(kind=8), dimension(:), allocatable, intent(inout) :: y
        character, optional :: trans

        character :: op
        integer :: uplo
        integer, dimension(:), allocatable :: rowind
        integer :: ii, jj, nnz


        if ( .not.present(trans) ) then
            op = 'N'
        else
            op = trans
        end if

        if ( op.eq.'N' .or. op.eq.'n' ) then
            if ( size(x).ne.A%cols ) then
                write(*,*) 'size(x) should be equal to A%cols in mso_sparse_z_mv!'
                stop 1
            end if
        else
            if ( op.eq.'c' .or. op.eq.'C' .or. op.eq.'t' .or. op.eq.'T' ) then
                if ( size(x).ne.A%rows ) then
                    write(*,*) 'size(x) should be equal to A%rows in mso_sparse_z_mv when trans is set to be N or n!'
                    stop 1
                end if
            else
                write(*,*) 'trans is not well defined in mso_sparse_z_mv!'
                stop 1
            end if
        end if

        if ( op.eq.'N' .or. op.eq.'n' ) then
            call cf_allocate(y,A%rows)
        else
            call cf_allocate(y,A%cols)
        end if

        y = 0.d0

        if ( op.eq.'N' .or. op.eq.'n' ) then
            do ii = 1, A%rows, 1
                do jj = A%ia(ii), A%ia(ii+1)-1, 1
                    y(ii) = y(ii) + A%a(jj)*x(A%ja(jj))
                end do
            end do
            return
        end if

        nnz = size(A%a)

        call cf_allocate(rowind,nnz)
        rowind = 0
        do ii = 1, A%rows, 1
            do jj = A%ia(ii), A%ia(ii+1)-1, 1
                rowind(jj) = ii
            end do
        end do

        if ( op.eq.'T' .or. op.eq.'t' ) then
            do ii = 1, nnz, 1
                y(A%ja(ii)) = y(A%ja(ii)) + A%a(ii)*x(rowind(ii))
            end do
        else
            do ii = 1, nnz, 1
                y(A%ja(ii)) = y(A%ja(ii)) + dconjg(A%a(ii))*x(rowind(ii))
            end do
        end if

        call cf_deallocate(rowind)
    end subroutine mso_sparse_z_mv

    subroutine mso_sparse_d_mm(A,B,C,trans)
        ! C = trans(A)*B
        ! trans = 'N', 'n', 'T', or 't'
        use common_functions, only: cf_allocate, cf_deallocate
        implicit none
        type(mso_dsparse_matrix), intent(in) :: A, B
        type(mso_dsparse_matrix), intent(inout) :: C
        character(len=1), intent(in), optional :: trans

        character(len=1) :: op
        integer :: nzmax, request, sort, m, n, k, info
        integer :: nnz

        if ( .not.present(trans) ) then
            op = 'N'
        else
            op = trans
        end if

        if ( op.eq.'N' .or. op.eq.'n' ) then
            if ( A%cols.ne.B%rows ) then
                write(*,*) 'The number of columns of A should be equal to the number of rows of B in mso_sparse_d_mm'
                stop 1
            end if
        else
            if ( op.eq.'T' .or. op.eq.'t' ) then
                if ( A%rows.ne.B%rows ) then
                    write(*,*) 'The number of rows of A should be equal to the number of rows of B in mso_sparse_d_mm'
                    stop 1
                end if
            else
                write(*,*) 'the input~trans is not well defined, which should be t, T, N, or n in mso_sparse_d_mm'
                stop 1
            end if
        end if

        m = A%rows
        n = A%cols
        k = B%rows
        sort = 9

        if ( op.eq.'T' .or. op.eq.'t' ) then
            C%rows = n
        else
            C%rows = m
        end if
        C%cols = B%cols

        call cf_allocate(C%ia,C%rows+1)

        request = 1
        call mkl_dcsrmultcsr(op, request, sort, m, n, k, A%a, A%ja, A%ia, B%a, B%ja, B%ia, C%a, C%ja, C%ia, nzmax, info)

        nnz = C%ia(C%rows+1) - 1

        if ( .not.allocated(C%a) .or. (size(C%a).lt.nnz) ) then
            call cf_allocate(C%a,nnz)
        end if
        if ( .not.allocated(C%ja) .or. (size(C%ja).lt.nnz) ) then
            call cf_allocate(C%ja,nnz)
        end if

        request = 2
        call mkl_dcsrmultcsr(op, request, sort, m, n, k, A%a, A%ja, A%ia, B%a, B%ja, B%ia, C%a, C%ja, C%ia, nzmax, info)

        if ( info.ne.0 ) then
            write(unit=*, fmt=*) 'Some wrong accurs in mso_sparse_d_mm: ifno:', info
            stop 1
        end if
    end subroutine mso_sparse_d_mm

    subroutine mso_sparse_d_mm2(A,B,trans)
        ! A = trans(A)*B
        implicit none
        type(mso_dsparse_matrix), intent(inout) :: A
        type(mso_dsparse_matrix), intent(in) :: B
        character(len=1), intent(in), optional :: trans

        type(mso_dsparse_matrix) :: C
        character(len=1) :: op

        if ( present(trans) ) then
            op = trans
        else
            op = 'N'
        end if

        call mso_sparse_d_mm(A,B,C,op)
        call mso_sparse_d_copy(C,A)

        call mso_sparse_destroy(C)
    end subroutine mso_sparse_d_mm2

    subroutine mso_sparse_z_mm(A,B,C,trans)
        ! C = trans(A)*B
        ! trans = 'N', 'n', 'T', or 't'
        use common_functions, only: cf_allocate
        type(mso_zsparse_matrix), intent(in) :: A, B
        type(mso_zsparse_matrix), intent(inout) :: C
        character(len=1), intent(in), optional :: trans

        character(len=1) :: op
        integer :: nzmax, request, sort, m, n, k, info
        integer :: nnz

        if ( .not.present(trans) ) then
            op = 'N'
        else
            op = trans
        end if

        if ( op.eq.'N' .or. op.eq.'n' ) then
            if ( A%cols.ne.B%rows ) then
                write(*,*) 'The number of columns of A should be equal to the number of rows of B in mso_sparse_z_mm'
                stop 1
            end if
        else
            if ( op.eq.'T' .or. op.eq.'t' .or. op.eq.'C' .or. op.eq.'c' ) then
                if ( A%rows.ne.B%rows ) then
                    write(*,*) 'The number of rows of A should be equal to the number of rows of B in mso_sparse_z_mm'
                    stop 1
                end if
            else
                write(*,*) 'the input~trans is not well defined, which should be t, T, C, c, N, or n in mso_sparse_z_mm'
                stop 1
            end if
        end if

        m = A%rows
        n = A%cols
        k = B%rows
        sort = 3

        if ( op.eq.'T' .or. op.eq.'t' .or. op.eq.'C' .or. op.eq.'c' ) then
            C%rows = n
        else
            C%rows = m
        end if
        C%cols = B%cols

        call cf_allocate(C%ia,C%rows+1)

        request = 1
        call mkl_zcsrmultcsr(op, request, sort, m, n, k, A%a, A%ja, A%ia, B%a, B%ja, B%ia, C%a, C%ja, C%ia, nzmax, info)

        nnz = C%ia(C%rows+1) - 1
        call cf_allocate(C%a,nnz)
        call cf_allocate(C%ja,nnz)

        request = 2
        call mkl_zcsrmultcsr(op, request, sort, m, n, k, A%a, A%ja, A%ia, B%a, B%ja, B%ia, C%a, C%ja, C%ia, nzmax, info)

        if ( info.ne.0 ) then
            write(unit=*, fmt=*) 'Some wrong accurs in mso_sparse_d_mm: ifno:', info
            stop 1
        end if
    end subroutine mso_sparse_z_mm

    subroutine mso_sparse_z_mm2(A,B,trans)
        ! A = trans(A)*B
        implicit none
        type(mso_zsparse_matrix), intent(inout) :: A
        type(mso_zsparse_matrix), intent(in) :: B
        character(len=1), intent(in), optional :: trans

        type(mso_zsparse_matrix) :: C
        character(len=1) :: op

        if ( present(trans) ) then
            op = trans
        else
            op = 'N'
        end if

        call mso_sparse_z_mm(A,B,C,op)
        call mso_sparse_z_copy(C,A)

        call mso_sparse_z_destroy(C)
    end subroutine mso_sparse_z_mm2

    subroutine mso_sparse_d_uppertriangle(upper,A,B)
        ! return the upper triangular part of A
        ! refered to triu in matlab
        use common_functions, only: cf_allocate, cf_deallocate
        implicit none
        integer, intent(in) :: upper
        type(mso_dsparse_matrix), intent(in) :: A
        type(mso_dsparse_matrix), intent(inout) :: B

        integer :: nnz, ii, jj, kk

        if ( .not.allocated(A%a) .or. .not.allocated(A%ia) .or. .not.allocated(A%ja) ) then
            return
        end if

        !         call mso_sparse_d_istriangle(A,ii)
        !         if ( ii.gt.0 ) then
        !             call mso_sparse_d_copy(A,B)
        !             return
        !         end if

        call cf_allocate(B%ia,A%rows+1)

        B%ia = 0
        do ii = 1, A%rows, 1
            do jj = A%ia(ii), A%ia(ii+1)-1, 1
                if ( A%ja(jj)-ii.ge.upper ) then
                    B%ia(ii+1) = B%ia(ii+1) + 1
                end if
            end do
        end do
        B%ia(1) = 1
        do ii = 2, A%rows+1, 1
            B%ia(ii) = B%ia(ii-1) + B%ia(ii)
        end do

        nnz = B%ia(A%rows+1) - 1
        if ( .not.allocated(B%a) .or. (size(B%a).lt.nnz) ) then
            call cf_allocate(B%a,nnz)
        end if
        if ( .not.allocated(B%ja) .or. (size(B%ja).lt.nnz) ) then
            call cf_allocate(B%ja,nnz)
        end if

        kk = 1
        do ii = 1, A%rows, 1
            do jj = A%ia(ii), A%ia(ii+1)-1, 1
                if ( A%ja(jj)-ii.ge.upper ) then
                    B%a(kk) = A%a(jj)
                    B%ja(kk) = A%ja(jj)
                    kk = kk + 1
                end if
            end do
        end do
        B%rows = A%rows
        B%cols = A%cols
    end subroutine mso_sparse_d_uppertriangle

    subroutine mso_sparse_z_uppertriangle(upper,A,B)
        ! return the upper triangular part of A
        ! refered to triu in matlab
        implicit none
        integer, intent(in) :: upper
        type(mso_zsparse_matrix), intent(in) :: A
        type(mso_zsparse_matrix), intent(inout) :: B

        integer :: nnz, ii, jj, kk

        if ( .not.allocated(A%a) .or. .not.allocated(A%ia) .or. .not.allocated(A%ja) ) then
            return
        end if

        call mso_sparse_z_istriangle(A,ii)
        if ( ii.gt.0 ) then
            call mso_sparse_z_copy(A,B)
            return
        end if

        if ( allocated(B%ia) ) then
            if ( size(B%ia).ne.A%rows+1 ) then
                deallocate(B%ia)
                allocate(B%ia(A%rows+1))
            end if
        else
            allocate(B%ia(A%rows+1))
        end if

        B%ia = 0
        do ii = 1, A%rows, 1
            do jj = A%ia(ii), A%ia(ii+1)-1, 1
                if ( A%ja(jj)-ii.ge.upper ) then
                    B%ia(ii+1) = B%ia(ii+1) + 1
                end if
            end do
        end do
        B%ia(1) = 1
        do ii = 2, A%rows+1, 1
            B%ia(ii) = B%ia(ii-1) + B%ia(ii)
        end do

        nnz = B%ia(A%rows+1) - 1
        if ( allocated(B%a) ) then
            if ( size(B%a).ne.nnz ) then
                deallocate(B%a)
                allocate(B%a(nnz))
            end if
        else
            allocate(B%a(nnz))
        end if
        if ( allocated(B%ja) ) then
            if ( size(B%ja).ne.nnz ) then
                deallocate(B%ja)
                allocate(B%ja(nnz))
            end if
        else
            allocate(B%ja(nnz))
        end if
        kk = 1
        do ii = 1, A%rows, 1
            do jj = A%ia(ii), A%ia(ii+1)-1, 1
                if ( A%ja(jj)-ii.ge.upper ) then
                    B%a(kk) = A%a(jj)
                    B%ja(kk) = A%ja(jj)
                    kk = kk + 1
                end if
            end do
        end do
        B%rows = A%rows
        B%cols = A%cols
    end subroutine mso_sparse_z_uppertriangle

    subroutine mso_sparse_d_istriangle(A,uplo)
        ! uplo: 1 if A is upper triangular
        !        -1 if A is lower triangular
        !        0 if A is not triangular
        type(mso_dsparse_matrix), intent(in) :: A
        integer, intent(inout) :: uplo

        integer :: ii, jj

        if ( .not.(allocated(A%a).and.allocated(A%ia).and.allocated(A%ja)) ) then
            write(*,*) 'A is a sparse matrix without non-zeros in mso_sparse_d_istriangle!'
            stop 1
        end if

        uplo = 1
        do ii = 1, A%rows, 1
            do jj = A%ja(A%ia(ii)), A%ja(A%ia(ii+1))-1, 1
                if ( A%ja(jj).lt.ii ) then
                    uplo = -1
                    exit
                end if
            end do
            if ( uplo.lt.0 ) exit
        end do
        if ( uplo.gt.0 ) return
        do ii = 1, A%rows, 1
            do jj = A%ja(A%ia(ii)), A%ja(A%ia(ii+1))-1, 1
                if ( A%ja(jj).gt.ii ) then
                    uplo = 0
                    exit
                end if
            end do
            if ( uplo.eq.0 ) exit
        end do
    end subroutine mso_sparse_d_istriangle

    subroutine mso_sparse_z_istriangle(A,uplo)
        ! uplo: 1 if A is upper triangular
        !        -1 if A is lower triangular
        !        0 if A is not triangular
        type(mso_zsparse_matrix), intent(in) :: A
        integer, intent(inout) :: uplo

        integer :: ii, jj

        if ( .not.(allocated(A%a).and.allocated(A%ia).and.allocated(A%ja)) ) then
            write(*,*) 'A is a sparse matrix without non-zeros in mso_sparse_d_istriangle!'
            stop 1
        end if

        uplo = 1
        do ii = 1, A%rows, 1
            do jj = A%ja(A%ia(ii)), A%ja(A%ia(ii+1))-1, 1
                if ( A%ja(jj).lt.ii ) then
                    uplo = -1
                    exit
                end if
            end do
            if ( uplo.lt.0 ) exit
        end do
        if ( uplo.gt.0 ) return
        do ii = 1, A%rows, 1
            do jj = A%ja(A%ia(ii)), A%ja(A%ia(ii+1))-1, 1
                if ( A%ja(jj).gt.ii ) then
                    uplo = 0
                    exit
                end if
            end do
            if ( uplo.eq.0 ) exit
        end do
    end subroutine mso_sparse_z_istriangle

    subroutine mso_sparse_d_pardiso(A,rhs,sol,msglvl)
        ! msglvl: 0 for no information output and 1 for information output
        use common_functions, only: cf_allocate
        implicit none
        type(mso_dsparse_matrix), intent(in) :: A
        real(kind=8), dimension(:), intent(inout) :: rhs
        real(kind=8), dimension(:), allocatable, intent(inout) :: sol
        integer, optional :: msglvl

        integer :: msg
        type(mso_dsparse_matrix) :: Smtx



        if ( A%rows.ne.A%cols ) then
            write(*,*) 'A%rows should be equal to A%cols in mso_sparse_d_pardiso'
            stop 1
        end if

        if ( A%rows.ne.size(rhs) ) then
            write(*,*) 'A%rows should be equal to size(rhs) in mso_sparse_d_pardiso'
            stop 1
        end if
        write(*,*) A%rows, allocated(sol), 'd pardiso start'
        call cf_allocate(sol,A%rows)
        write(*,*) 'd pardiso xxxx'
        if ( present(msglvl) ) then
            msg = msglvl ! msglvl: 0 for no information output and 1 for information output
        else
            msg = 0
        end if

        if ( A%mtype.eq.2 .or. A%mtype.eq.-2 ) then
            call mso_sparse_uppertriangle(0,A,Smtx)
            call mso_dpardiso(A%mtype,msg,Smtx%a,Smtx%ia,Smtx%ja,rhs,sol)
        else
            call mso_dpardiso(A%mtype,msg,A%a,A%ia,A%ja,rhs,sol)
        end if

        call mso_sparse_destroy(Smtx)
    end subroutine mso_sparse_d_pardiso

    subroutine mso_sparse_z_pardiso(A,rhs,x,msglvl)
        ! msglvl: 0 for no information output and 1 for information output
        implicit none
        type(mso_zsparse_matrix), intent(in) :: A
        complex(kind=8), dimension(:), intent(inout) :: rhs
        complex(kind=8), dimension(:), allocatable, intent(inout) :: x
        integer, optional :: msglvl

        integer :: msg
        type(mso_zsparse_matrix) :: Smtx

        if ( A%rows.ne.A%cols ) then
            write(*,*) 'A%rows should be equal to A%cols in mso_sparse_d_pardiso'
            stop 1
        end if

        if ( A%rows.ne.size(rhs) ) then
            write(*,*) 'A%rows should be equal to size(rhs) in mso_sparse_d_pardiso'
            stop 1
        end if

        if ( allocated(x) ) then
            if ( size(x).ne.A%cols ) then
                deallocate(x)
                allocate(x(A%cols))
            end if
        else
            allocate(x(A%cols))
        end if

        if ( present(msglvl) ) then
            msg = msglvl ! msglvl: 0 for no information output and 1 for information output
        else
            msg = 0
        end if

        if ( A%mtype.eq.4 .or. A%mtype.eq.-4 ) then
            call mso_sparse_uppertriangle(0,A,Smtx)
            call mso_zpardiso(A%mtype,msg,Smtx%a,Smtx%ia,Smtx%ja,rhs,x)
        else
            call mso_zpardiso(A%mtype,msg,A%a,A%ia,A%ja,rhs,x)
        end if

        call mso_sparse_destroy(Smtx)
    end subroutine mso_sparse_z_pardiso

    subroutine mso_sparse_d_gmres(maxiter, maxrelres, A, b, x, itercount, residuals,msglvl)
        !*******************************************************************************
        !  Content:
        !  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
        !  to solve the linear system A*x = b to get the solution x
        !  msglvl: print the procesure data if msglvl>0
        !  residuals(2): residuals(1) is the Euclidean norm of the current residual
        !                residuals(2) is equal to res(1) / norm2(rhs)
        !*******************************************************************************
        implicit none
        integer, intent(in) :: maxiter(2)
        real(kind=8), intent(in) :: maxrelres
        type(mso_dsparse_matrix), intent(in) :: A
        real(kind=8), dimension(:), intent(in) :: b
        real(kind=8), dimension(:), allocatable, intent(inout) :: x
        integer, intent(inout) :: itercount
        real(kind=8), intent(inout) :: residuals(2)
        integer, optional :: msglvl

        integer :: msg

        if ( present(msglvl) ) then
            msg = msglvl
        else
            msg = -1
        end if

        if ( A%rows.ne.A%cols ) then
            write(*,*) 'A is not a square matrix in mso_sparse_d_gmres!'
            stop 1
        end if

        if ( A%rows.ne.size(b) ) then
            write(*,*) 'A%rows is not equal to size(b) in mso_sparse_d_gmres!'
            stop 1
        end if

        if ( allocated(x) ) then
            if ( size(x).ne.A%cols ) then
                deallocate(x)
                allocate(x(A%cols))
                x = 0.d0 ! set initial solution
            end if
        else
            allocate(x(A%cols))
            x = 0.d0 ! set initial solution
        end if

        call mso_dgmres_nonprecondition(maxiter, maxrelres, A%a,A%ia,A%ja,b, x, itercount, residuals,msg)
    end subroutine mso_sparse_d_gmres


    subroutine mso_sparse_d_show(A, fmt, trans, isfull)
        ! trans: 'N', 'n', 'T', 't', 'C', 'c'
        ! isfull: not pesent or >0 if full, <0 otherwise
        use common_functions, only: cf_allocate, cf_deallocate
        implicit none
        type(mso_dsparse_matrix), intent(in) :: A
        character(len=*), intent(in) :: fmt
        character(len=1), intent(in), optional :: trans
        integer, intent(in), optional :: isfull

        integer :: ii, jj, kk
        character(len=20) :: s1
        type(mso_dsparse_matrix) :: B

        if ( A%rows*A%cols.eq.0 ) return

        if ( .not.present(isfull) .or. isfull.gt.0 ) then
            if ( .not.present(trans) .or. trans.eq.'N' .or. trans.eq.'n' ) then
                do ii = 1, A%rows, 1
                    if ( A%ia(ii+1)-A%ia(ii).eq.0 ) then
                        do jj = 1, A%cols, 1
                            write(*,'('//fmt//'$)') 0.d0
                        end do
                    else
                        kk = A%ia(ii)
                        do jj = 1, A%cols, 1
                            if ( jj.eq.A%ja(kk) ) then
                                write(*,'('//fmt//'$)') A%a(kk)
                                if ( kk.lt.A%ia(ii+1)-1 ) kk = kk + 1
                            else
                                write(*,'('//fmt//'$)') 0.d0
                            end if
                        end do
                    end if
                    write(*,*)
                end do
                return
            end if

            if ( trans.eq.'T' .or. trans.eq.'t' .or. trans.eq.'c' .or. trans.eq.'C' ) then
                call mso_sparse_d_transpose(A,B,'t')
                do ii = 1, B%rows, 1
                    if ( B%ia(ii+1)-B%ia(ii).eq.0 ) then
                        do jj = 1, B%cols, 1
                            write(*,'('//fmt//'$)') 0.d0
                        end do
                    else
                        kk = B%ia(ii)
                        do jj = 1, B%cols, 1
                            if ( jj.eq.B%ja(kk) ) then
                                write(*,'('//fmt//'$)') B%a(kk)
                                if ( kk.lt.B%ia(ii+1)-1 ) kk = kk + 1
                            else
                                write(*,'('//fmt//'$)') 0.d0
                            end if
                        end do
                    end if
                    write(*,*)
                end do
            else
                write(*,*) 'trans is not well defined in mso_sparse_d_show!'
                stop 1
            end if
        else
            if ( .not.present(trans) .or. trans.eq.'N' .or. trans.eq.'n' ) then
                do ii = 1, A%rows, 1
                    do jj = A%ia(ii), A%ia(ii+1)-1, 1
                        write(*,'(I5,I5,'//fmt//')') ii, A%ja(jj), A%a(jj)
                    end do
                end do
            else
                do ii = 1, A%rows, 1
                    do jj = A%ia(ii), A%ia(ii+1)-1, 1
                        write(*,'(I5,I5,'//fmt//')') A%ja(jj), ii, A%a(jj)
                    end do
                end do
            end if
        end if

        call mso_sparse_d_destroy(B)
    end subroutine mso_sparse_d_show

    subroutine mso_sparse_z_show(A, fmt, isfull)
        type(mso_zsparse_matrix), intent(in) :: A
        character(len=*), intent(in) :: fmt
        integer, intent(in), optional :: isfull

        integer :: ii, jj
        character(len=20) :: s1, s2

        if ( A%rows*A%cols.eq.0 ) then
            write(*,*) 'No nonzero elements exist in the double complex sparse matrix!'
            return
        end if

        if (present(isfull) .and. isfull.gt.0 ) then
            call mso_sparse_showZ(A%a, A%ia, A%ja, (/1,A%rows/), (/1,A%cols/), fmt)
            return
        end if
        do ii = 1, A%rows, 1
            do jj = A%ia(ii), A%ia(ii+1)-1, 1
                write(s1,fmt) dreal(A%a(jj))
                if ( dimag(A%a(jj)).ge.0.d0 ) then
                    write(s2,fmt) dimag(A%a(jj))
                    write(*,*) ii, A%ja(jj), '    '//trim(adjustl(s1))//'+'//trim(adjustl(s2))//'i'
                else
                    write(s2,fmt) -dimag(A%a(jj))
                    write(*,*) ii, A%ja(jj), '    '//trim(adjustl(s1))//'-'//trim(adjustl(s2))//'i'
                end if
            end do
        end do
    end subroutine mso_sparse_z_show



    subroutine mso_sparse_d_destroy(A)
        type(mso_dsparse_matrix), intent(inout) :: A

        integer :: err

        if (allocated(A%a)) then
            deallocate(A%a, stat=err)
            if (err /= 0) print *, "A%a in mso_allocateDSpMtx: Deallocation request denied"
        end if
        if (allocated(A%ia)) then
            deallocate(A%ia, stat=err)
            if (err /= 0) print *, "A%ia in mso_allocateDSpMtx: Deallocation request denied"
        end if
        if (allocated(A%ja)) then
            deallocate(A%ja, stat=err)
            if (err /= 0) print *, "A%ja in mso_allocateDSpMtx: Deallocation request denied"
        end if
        A%rows = 0
        A%cols = 0
    end subroutine mso_sparse_d_destroy

    subroutine mso_sparse_z_destroy(A)
        type(mso_zsparse_matrix), intent(inout) :: A

        integer :: err

        if (allocated(A%a)) then
            deallocate(A%a, stat=err)
            if (err /= 0) print *, "A%a in mso_allocateDSpMtx: Deallocation request denied"
        end if
        if (allocated(A%ia)) then
            deallocate(A%ia, stat=err)
            if (err /= 0) print *, "A%ia in mso_allocateDSpMtx: Deallocation request denied"
        end if
        if (allocated(A%ja)) then
            deallocate(A%ja, stat=err)
            if (err /= 0) print *, "A%ja in mso_allocateDSpMtx: Deallocation request denied"
        end if
        A%rows = 0
        A%cols = 0
    end subroutine mso_sparse_z_destroy



    ! old codes



    subroutine mso_zgmres_precondition(maxiter, maxrelres, a,ia,ja,rhs, computed_solution, itercount, relres, prefun)
        !*******************************************************************************
        !  Content:
        !  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
        !*******************************************************************************
        implicit none
        interface
            subroutine prefun(x,u)
                complex(kind=8), dimension(:), intent(in) :: x
                complex(kind=8), dimension(:), intent(inout) :: u
            end subroutine prefun
        end interface
        integer, intent(in) :: maxiter(2)
        real(kind=8), intent(in) :: maxrelres
        complex(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        complex(kind=8), dimension(:), intent(in) :: rhs
        complex(kind=8), dimension(:), intent(inout) :: computed_solution
        integer, intent(inout) :: itercount
        real(kind=8), intent(inout) :: relres


        integer :: n, zn
        integer :: ipar(128)
        real(kind=8) :: dpar(128)
        real(kind=8), dimension(:), allocatable :: tmp, b, drhs, dcs
        complex(kind=8), dimension(:), allocatable :: residual

        integer :: rci_request, i
        real(kind=8) :: dvar, rhsm2
        real(kind=8), external :: dznrm2
        integer :: err

        err = 0

        zn = size(ia) - 1
        n = 2*zn    ! the size of the problem
        rhsm2 = dznrm2(zn,rhs,1)

        allocate( b(n))
        allocate( residual(zn))
        allocate( drhs(n))
        allocate( dcs(n))

        !---------------------------------------------------------------------------
        ! Initialize variables and the right hand side through matrix-vector product
        !---------------------------------------------------------------------------
        call dcopy(zn, dreal(rhs), 1, drhs, 1)
        call dcopy(zn, dimag(rhs), 1, drhs(zn+1), 1)

        !       dcs = 0.d0
        call prefun(rhs,computed_solution)
        call dcopy(zn, dreal(computed_solution), 1, dcs, 1)
        call dcopy(zn, dimag(computed_solution), 1, dcs(zn+1), 1)

        ipar(15) = maxiter(2)
        allocate( tmp(n*(2*ipar(15)+1)+(ipar(15)*(ipar(15)+9))/2+1))

        call dfgmres_init(n, dcs, drhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        ipar(8) = 1
        ipar(15) = maxiter(2)
        ipar(5) = maxiter(1)
        ipar(11) = 1    ! precondition
        dpar(1) = maxrelres

        call dfgmres_check(n, dcs, drhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        1 call dfgmres(n, dcs, drhs, rci_request, ipar, dpar, tmp)

        if (rci_request .eq.0) go to 3

        if (rci_request .eq. 1) then
            computed_solution = dcmplx( tmp(ipar(22):ipar(22)+zn-1),    &
                    &   tmp(ipar(22)+zn:ipar(23)-1) )
        !           call zcopy(zn,dcmplx(tmp(ipar(22):ipar(22)+zn-1),tmp(ipar(22)+zn:ipar(23)-1)),1,computed_solution,1)
            call mkl_zcsrgemv('n',zn, a, ia, ja, computed_solution, residual)
            call dcopy(zn,dreal(residual),1,tmp(ipar(23)),1)
            call dcopy(zn,dimag(residual),1, tmp(ipar(23)+zn), 1)
            go to 1
        end if

        if (rci_request .eq. 2) then
            ! request to the dfgmres_get routine to put the solution into b(n) via ipar(13)
            ipar(13) = 1
            ! get the current fgmres solution in the vector b(n)
            call dfgmres_get(n, dcs, b, rci_request, ipar,dpar, tmp, itercount)
            ! compute the current true residual via mkl (sparse) blas routines
            computed_solution = dcmplx(b(1:zn),b(zn+1:n))
            call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
            call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
            dvar = dznrm2(zn, residual, 1)
            dvar = dvar/rhsm2
            if (dvar .lt. maxrelres) then
                go to 3
            else
                go to 1
            end if
        end if

        if (rci_request .eq. 3) then
            computed_solution = dcmplx( tmp(ipar(22):ipar(22)+zn-1),    &
                    &   tmp(ipar(22)+zn:ipar(23)-1) )
            call prefun( computed_solution, residual )
            call dcopy(zn, dreal(residual), 1, tmp(ipar(23)), 1)
            call dcopy(zn, dimag(residual), 1, tmp(ipar(23)+zn), 1)
            go to 1
        end if

        if (rci_request .eq. 4) then
            if (dpar(7) .lt. 1.0d-12) then
                go to 3
            else
                go to 1
            end if
        else
            go to 999
        end if

        3 ipar(13) = 0
        call dfgmres_get(n, dcs, drhs, rci_request, ipar, dpar, tmp, itercount)
        computed_solution = dcmplx( dcs(1:zn), dcs(zn+1:n) )
        call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
        call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
        relres = dznrm2(zn, residual, 1)
        relres = relres/rhsm2

        if (allocated(dcs)) deallocate(dcs, stat=err)
        if (err /= 0) print *, "dcs: Deallocation request denied"
        if (allocated(drhs)) deallocate(drhs, stat=err)
        if (err /= 0) print *, "drhs: Deallocation request denied"
        if (allocated(tmp)) deallocate(tmp, stat=err)
        if (err /= 0) print *, "tmp: Deallocation request denied"
        if (allocated(b)) deallocate(b, stat=err)
        if (err /= 0) print *, "b: Deallocation request denied"
        if (allocated(residual)) deallocate(residual, stat=err)
        if (err /= 0) print *, "residual: Deallocation request denied"
        return


        999 write( *,'(A,I5)') 'this example failed as the solver has returned the error code', rci_request
        call mkl_free_buffers
        stop 1
    end subroutine mso_zgmres_precondition

    subroutine mso_zgmres_nonprecondition(maxiter, maxrelres, a,ia,ja,rhs, computed_solution, itercount, relres)
        !*******************************************************************************
        !  Content:
        !  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
        !*******************************************************************************
        implicit none
        integer, intent(in) :: maxiter(2)
        real(kind=8), intent(in) :: maxrelres
        complex(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        complex(kind=8), dimension(:), intent(in) :: rhs
        complex(kind=8), dimension(:), intent(inout) :: computed_solution
        integer, intent(inout) :: itercount
        real(kind=8), intent(inout) :: relres


        integer :: n, zn
        integer :: ipar(128)
        real(kind=8) :: dpar(128)
        real(kind=8), dimension(:), allocatable :: tmp, b, drhs, dcs
        complex(kind=8), dimension(:), allocatable :: residual

        integer :: rci_request, i
        real(kind=8) :: dvar, rhsm2
        real(kind=8), external :: dznrm2
        integer :: err

        err = 0

        zn = size(ia) - 1
        n = 2*zn    ! the size of the problem
        rhsm2 = dznrm2(zn,rhs,1)

        allocate( b(n))
        allocate( residual(zn))
        allocate( drhs(n))
        allocate( dcs(n))

        !---------------------------------------------------------------------------
        ! Initialize variables and the right hand side through matrix-vector product
        !---------------------------------------------------------------------------
        call dcopy(zn, dreal(rhs), 1, drhs, 1)
        call dcopy(zn, dimag(rhs), 1, drhs(zn+1), 1)

        dcs = 1.d0

        ipar(15) = maxiter(2)
        allocate( tmp(n*(2*ipar(15)+1)+(ipar(15)*(ipar(15)+9))/2+1))

        call dfgmres_init(n, dcs, drhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        ipar(8) = 1
        ipar(15) = maxiter(2)
        ipar(5) = maxiter(1)
        ipar(11) = 0
        dpar(1) = maxrelres

        call dfgmres_check(n, dcs, drhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        1 call dfgmres(n, dcs, drhs, rci_request, ipar, dpar, tmp)

        if (rci_request .eq.0) go to 3

        if (rci_request .eq. 1) then
            computed_solution = dcmplx( tmp(ipar(22):ipar(22)+zn-1),    &
                    &   tmp(ipar(22)+zn:ipar(23)-1) )
            call mkl_zcsrgemv('n',zn, a, ia, ja, computed_solution, residual)
            call dcopy(zn,dreal(residual),1,tmp(ipar(23)),1)
            call dcopy(zn,dimag(residual),1, tmp(ipar(23)+zn), 1)
            go to 1
        end if

        if (rci_request .eq. 2) then
            ! request to the dfgmres_get routine to put the solution into b(n) via ipar(13)
            ipar(13) = 1
            ! get the current fgmres solution in the vector b(n)
            call dfgmres_get(n, dcs, b, rci_request, ipar,dpar, tmp, itercount)
            ! compute the current true residual via mkl (sparse) blas routines
            computed_solution = dcmplx(b(1:zn),b(zn+1:n))
            call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
            call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
            dvar = dznrm2(zn, residual, 1)
            dvar = dvar/rhsm2
            if (dvar .lt. maxrelres) then
                go to 3
            else
                go to 1
            end if
        end if

        if (rci_request .eq. 3) then
            write(unit=*, fmt=*) 3

            call dcopy(n,tmp(ipar(22)),1, tmp(ipar(23)), 1)
            go to 1
        end if

        if (rci_request .eq. 4) then
            if (dpar(7) .lt. 1.0d-12) then
                go to 3
            else
                go to 1
            end if
        else
            go to 999
        end if

        3 ipar(13) = 0
        call dfgmres_get(n, dcs, drhs, rci_request, ipar, dpar, tmp, itercount)
        computed_solution = dcmplx( dcs(1:zn), dcs(zn+1:n) )
        call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
        call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
        relres = dznrm2(zn, residual, 1)
        relres = relres/rhsm2

        if (allocated(dcs)) deallocate(dcs, stat=err)
        if (err /= 0) print *, "dcs: Deallocation request denied"
        if (allocated(drhs)) deallocate(drhs, stat=err)
        if (err /= 0) print *, "drhs: Deallocation request denied"
        if (allocated(tmp)) deallocate(tmp, stat=err)
        if (err /= 0) print *, "tmp: Deallocation request denied"
        if (allocated(b)) deallocate(b, stat=err)
        if (err /= 0) print *, "b: Deallocation request denied"
        if (allocated(residual)) deallocate(residual, stat=err)
        if (err /= 0) print *, "residual: Deallocation request denied"
        return


        999 write( *,'(A,I5)') 'this example failed as the solver has returned the error code', rci_request
        call mkl_free_buffers
        stop 1
    end subroutine mso_zgmres_nonprecondition

    subroutine mso_dcoo2csr(rowind,colind,acoo,acsr,ia,ja,m,nnz)
        integer, intent(in) :: m,nnz
        integer, intent(in) :: rowind(nnz), colind(nnz)
        real(kind=8), intent(in) :: acoo(nnz)
        real(kind=8), intent(inout) :: acsr(nnz)
        integer, intent(inout) :: ia(m+1), ja(nnz)

        integer :: job(8), info

        job(1) = 2  ! converted to CSR
        job(2) = 1  ! one-based indexing for CSR
        job(3) = 1  ! one-based indexing for COO

        job(6) = 1
        call mkl_dcsrcoo( job, m, acsr, ja, ia, nnz, acoo, rowind, colind, info)
        job(6) = 2
        call mkl_dcsrcoo( job, m, acsr, ja, ia, nnz, acoo, rowind, colind, info)
    end subroutine mso_dcoo2csr

    subroutine mso_dcsradd(trans,m,n, acsr,ja,ia, beta, bcsr,jb,ib, ccsr,jc,ic)
        character, intent(in) :: trans
        integer, intent(in) :: m,n, ia(*), ja(*), ib(*), jb(*)
        real(kind=8), intent(in) :: beta, acsr(*), bcsr(*)
        real(kind=8), dimension(:), intent(inout) :: ccsr
        integer, dimension(:), intent(inout) :: jc, ic

        integer :: info, nzmax

        call mkl_dcsradd(trans,1,3,m,n, acsr,ja,ia, beta, bcsr,jb,ib, ccsr,jc,ic,nzmax, info)
        call mkl_dcsradd(trans,2,3,m,n, acsr,ja,ia, beta, bcsr,jb,ib, ccsr,jc,ic,nzmax, info)
    end subroutine mso_dcsradd

    subroutine mso_dcsrdltRC(a,ia,ja,v,RC, b,ib,jb, deletedcols)
        ! (a,ia,ja) CSR formart sparse matrix,
        ! v is a integer vector with elements > 0,
        ! (b,ib,jb) should not be allocated,
        ! Rc='R' or 'C'. 'R' represents row; 'C' represents column
        use common_functions, only: cf_qsort_i, cf_allocate, cf_deallocate
        implicit none
        real(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja, v
        character, intent(in) :: RC
        real(kind=8), dimension(:), allocatable, intent(inout) :: b
        integer, dimension(:), allocatable, intent(inout) :: ib,jb
        integer, optional :: deletedcols

        integer, dimension(:), allocatable :: vv, itmp
        real(kind=8), dimension(:), allocatable :: dtmp
        integer :: ii, m, n, sv, jj, sib, kk, nnz, err

        sv = size(v)
        allocate(vv(sv))
        vv = v
        call cf_qsort_i(vv,sv,mso_cmp)   ! sort vv
        ! delete repetitive numbers
        if ( sv.gt.1 ) then
            ii = 2
            do while ( ii.le.sv .and. vv(ii-1).ne.vv(ii) )
                ii = ii + 1
            end do
            jj = ii + 1
            do while ( jj.le.sv )
                do while ( jj.le.sv .and. vv(jj-1).eq.vv(jj) )
                    jj = jj + 1
                end do
                if ( jj.le.sv ) then
                    vv(ii) = vv(jj)
                    ii = ii + 1
                    jj = jj + 1
                end if
            end do
            sv = ii - 1
        end if

        if ( RC.eq.'R' .or. RC.eq.'r' ) then
            ! delete row numbers>m
            m = size(ia) - 1
            ii = sv
            do while ( ii.ge.1 .and. vv(ii).gt.m )
                ii = ii - 1
            end do
            sv = ii
            if ( present(deletedcols) ) then
                deletedcols = 0
            end if
            if ( sv.eq.0 ) then
                jj = size(ia)
                kk = ia(jj) - 1
                call cf_allocate(ib,jj)
                ib = ia
                if ( .not.allocated(jb) .or. (size(ib).lt.kk) ) then
                    call cf_allocate(jb,kk)
                end if
                ib(1:kk) = ia(1:kk)
                if ( .not.allocated(b) .or. (size(b).lt.kk) ) then
                    call cf_allocate(b,kk)
                end if
                b(1:kk) = a(1:kk)
                goto 100
            end if
            sib = m - sv + 1    ! size(ib)
            call cf_allocate(ib,sib)
            allocate(itmp(sib-1))
            jj = 1
            kk = 1
            nnz = 0
            ! value b, ib, jb
            do ii = 1, m, 1
                if ( jj.le.sv .and. ii.eq.vv(jj) ) then
                    jj = jj + 1
                else
                    itmp(kk) = ii
                    ib(kk) = nnz + 1
                    nnz = nnz + ia(ii+1) - ia(ii)
                    kk = kk + 1
                end if
            end do
            ib(kk) = nnz + 1
            if ( .not.allocated(jb) .or. (size(jb).lt.nnz) ) then
                call cf_allocate(jb,nnz)
            end if
            if ( .not.allocated(b) .or. (size(b).lt.nnz) ) then
                call cf_allocate(b,nnz)
            end if
            forall ( ii = 1:sib-1 )
                jb(ib(ii):ib(ii+1)-1) = ja(ia(itmp(ii)):ia(itmp(ii)+1)-1)
                b(ib(ii):ib(ii+1)-1) = a(ia(itmp(ii)):ia(itmp(ii)+1)-1)
            end forall
        else
            nnz = ia(size(ia)) - 1
            m = size(ia) - 1
            n = maxval(ja(1:nnz))
            ii = sv
            do while ( ii.ge.1 .and. vv(ii).gt.n )
                ii = ii - 1
            end do
            sv = ii
            if ( present(deletedcols) ) then
                deletedcols = sv
            end if
            if ( sv.eq.0 ) then
                jj = size(ia)
                kk = ia(jj) - 1
                call cf_allocate(ib,jj)
                ib = ia
                if ( .not.allocated(jb) .or. (size(ib).lt.kk) ) then
                    call cf_allocate(jb,kk)
                end if
                ib(1:kk) = ia(1:kk)
                if ( .not.allocated(b) .or. (size(b).lt.kk) ) then
                    call cf_allocate(b,kk)
                end if
                b(1:kk) = a(1:kk)
                goto 100
            end if
            call cf_allocate(ib, size(ia))
            call cf_allocate(itmp,nnz)
            call cf_allocate(dtmp,nnz)
            itmp = ja(1:nnz)
            ib = 0
            do ii = 1, sv, 1
                do jj = 1, m, 1
                    do kk = ia(jj), ia(jj+1)-1, 1
                        if ( ja(kk).gt.vv(ii) ) then
                            itmp(kk) = itmp(kk) - 1
                        else
                            if ( ja(kk).eq.vv(ii) ) then
                                itmp(kk) = 0
                                ib(jj+1) = ib(jj+1) + 1
                            end if
                        end if
                    end do
                end do
            end do
            ib(1) = 1
            do ii = 2, m+1, 1
                ib(ii) = ib(ii-1) + ia(ii) - ia(ii-1) - ib(ii)
            end do

            dtmp = a(1:nnz)
            ii = 1
            do while ( ii.le.nnz .and. itmp(ii).ne.0 )
                ii = ii + 1
            end do
            if ( ii.le.nnz ) then
                jj = ii + 1
                do while ( jj.le.nnz )
                    if ( itmp(jj).eq.0 ) then
                        jj = jj + 1
                    else
                        itmp(ii) = itmp(jj)
                        dtmp(ii) = dtmp(jj)
                        ii = ii + 1
                        jj = jj + 1
                    end if
                end do
            end if
            nnz = ii - 1
            if ( .not.allocated(b) .or. (size(b).lt.nnz) ) then
                call cf_allocate(b,nnz)
            end if
            b(1:nnz) = dtmp(1:nnz)
            if ( .not.allocated(jb) .or. (size(jb).lt.nnz) ) then
                call cf_allocate(jb,nnz)
            end if
            jb(1:nnz) = itmp(1:nnz)
        end if

        100 continue

        call cf_deallocate(itmp)
        call cf_deallocate(vv)
        call cf_deallocate(dtmp)
    contains
        function mso_cmp(a,b)
            integer, intent(in) :: a,b
            integer :: mso_cmp
            if ( a.lt.b ) mso_cmp = -1
            if ( a.eq.b ) mso_cmp = 0
            if ( a.gt.b ) mso_cmp = 1
        end function mso_cmp
    end subroutine mso_dcsrdltRC

    subroutine mso_zcsrdltRC(a,ia,ja,v,RC, b,ib,jb, deletedcols)
        ! (a,ia,ja) CSR formart sparse matrix,
        ! v is a integer vector with elements > 0,
        ! (b,ib,jb) should not be allocated,
        ! Rc='R' or 'C'. 'R' represents row; 'C' represents column
        use common_functions, only: cf_qsort_i
        implicit none
        complex(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja, v
        character, intent(in) :: RC
        complex(kind=8), dimension(:), allocatable :: b
        integer, dimension(:), allocatable :: ib,jb
        integer, optional :: deletedcols

        integer, dimension(:), allocatable :: vv, itmp
        complex(kind=8), dimension(:), allocatable :: dtmp
        integer :: ii, m, n, sv, jj, sib, kk, nnz, err

        if (allocated(b)) then
            deallocate(b, stat=err)
            if (err /= 0) print *, "b in mso_dcsrdltRC: Deallocation request denied"
        end if
        if (allocated(ib)) then
            deallocate(ib, stat=err)
            if (err /= 0) print *, "ib in mso_dcsrdltRC: Deallocation request denied"
        end if
        if (allocated(jb)) then
            deallocate(jb, stat=err)
            if (err /= 0) print *, "jb in mso_dcsrdltRC: Deallocation request denied"
        end if

        sv = size(v)
        allocate(vv(sv))
        vv = v
        call cf_qsort_i(vv,sv,mso_cmp)   ! sort vv
        ! delete repetitive numbers
        if ( sv.gt.1 ) then
            ii = 2
            do while ( ii.le.sv .and. vv(ii-1).ne.vv(ii) )
                ii = ii + 1
            end do
            jj = ii + 1
            do while ( jj.le.sv )
                do while ( jj.le.sv .and. vv(jj-1).eq.vv(jj) )
                    jj = jj + 1
                end do
                if ( jj.le.sv ) then
                    vv(ii) = vv(jj)
                    ii = ii + 1
                    jj = jj + 1
                end if
            end do
            sv = ii - 1
        end if

        if ( RC.eq.'R' .or. RC.eq.'r' ) then
            ! delete row numbers>m
            m = size(ia) - 1
            ii = sv
            do while ( ii.ge.1 .and. vv(ii).gt.m )
                ii = ii - 1
            end do
            sv = ii
            if ( present(deletedcols) ) then
                deletedcols = 0
            end if
            if ( sv.eq.0 ) then
                allocate(b(size(a)))
                b = a
                allocate(ib(size(ia)))
                ib = ia
                allocate(jb(size(ja)))
                jb = ja
                goto 100
            end if
            sib = m - sv + 1    ! size(ib)
            allocate(ib(sib))
            allocate(itmp(sib-1))
            jj = 1
            kk = 1
            nnz = 0
            ! value b, ib, jb
            do ii = 1, m, 1
                if ( jj.le.sv .and. ii.eq.vv(jj) ) then
                    jj = jj + 1
                else
                    itmp(kk) = ii
                    ib(kk) = nnz + 1
                    nnz = nnz + ia(ii+1) - ia(ii)
                    kk = kk + 1
                end if
            end do
            ib(kk) = nnz + 1
            allocate(jb(nnz))
            allocate(b(nnz))
            forall ( ii = 1:sib-1 )
                jb(ib(ii):ib(ii+1)-1) = ja(ia(itmp(ii)):ia(itmp(ii)+1)-1)
                b(ib(ii):ib(ii+1)-1) = a(ia(itmp(ii)):ia(itmp(ii)+1)-1)
            end forall
        else
            m = size(ia) - 1
            n = maxval(ja)
            ii = sv
            do while ( ii.ge.1 .and. vv(ii).gt.n )
                ii = ii - 1
            end do
            sv = ii
            if ( present(deletedcols) ) then
                deletedcols = sv
            end if
            if ( sv.eq.0 ) then
                allocate(b(size(a)))
                b = a
                allocate(ib(size(ia)))
                ib = ia
                allocate(jb(size(ja)))
                jb = ja
                goto 100
            end if
            allocate(ib(size(ia)))
            nnz = size(ja)
            if (allocated(itmp)) then
                deallocate(itmp, stat=err)
                if (err /= 0) print *, "itmp in mso_dcsrdltRC: Deallocation request denied"
            end if
            allocate(itmp(nnz), stat=err)
            if (err /= 0) print *, "itmp in mso_dcsrdltRC: Allocation request denied"
            allocate(dtmp(nnz), stat=err)
            if (err /= 0) print *, "dtmp in mso_dcsrdltRC: Allocation request denied"
            itmp = ja
            ib = 0
            do ii = 1, sv, 1
                do jj = 1, m, 1
                    do kk = ia(jj), ia(jj+1)-1, 1
                        if ( ja(kk).gt.vv(ii) ) then
                            itmp(kk) = itmp(kk) - 1
                        else
                            if ( ja(kk).eq.vv(ii) ) then
                                itmp(kk) = 0
                                ib(jj+1) = ib(jj+1) + 1
                            end if
                        end if
                    end do
                end do
            end do
            ib(1) = 1
            do ii = 2, m+1, 1
                ib(ii) = ib(ii-1) + ia(ii) - ia(ii-1) - ib(ii)
            end do

            dtmp = a
            ii = 1
            do while ( ii.le.nnz .and. itmp(ii).ne.0 )
                ii = ii + 1
            end do
            if ( ii.le.nnz ) then
                jj = ii + 1
                do while ( jj.le.nnz )
                    if ( itmp(jj).eq.0 ) then
                        jj = jj + 1
                    else
                        itmp(ii) = itmp(jj)
                        dtmp(ii) = dtmp(jj)
                        ii = ii + 1
                        jj = jj + 1
                    end if
                end do
            end if
            nnz = ii - 1
            allocate(b(nnz), stat=err)
            if (err /= 0) print *, "b in mso_dcsrdltRC: Allocation request denied"
            b = dtmp(1:nnz)
            allocate(jb(nnz), stat=err)
            if (err /= 0) print *, "jb in mso_dcsrdltRC: Allocation request denied"
            jb = itmp(1:nnz)
        end if

        100 continue

        if ( allocated(itmp) ) deallocate(itmp)
        if ( allocated(vv) ) deallocate(vv)
        if (allocated(dtmp)) then
            deallocate(dtmp, stat=err)
            if (err /= 0) print *, "dtmp in mso_dcsrdltRC: Deallocation request denied"
        end if
    contains
        function mso_cmp(a,b)
            integer, intent(in) :: a,b
            integer :: mso_cmp
            if ( a.lt.b ) mso_cmp = -1
            if ( a.eq.b ) mso_cmp = 0
            if ( a.gt.b ) mso_cmp = 1
        end function mso_cmp
    end subroutine mso_zcsrdltRC

    subroutine mso_sparse_showz(a, ia, ja, rows, cols, fmt, trans)
        use common_functions, only: prtVec
        implicit none
        complex(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        integer, intent(in) :: rows(2), cols(2)
        character(len=*), intent(in), optional :: fmt
        character(len=1), intent(in), optional :: trans

        integer :: i, mn(2), m, job(6), info, n
        complex(kind=8), dimension(:,:), allocatable :: adns

        m = size(ia) - 1
        n = maxval(ja)
        allocate(adns(m,n))

        job = 1
        job(4) = 2
        call mkl_zdnscsr(job,m,n,adns,m,a,ja,ia,info)

        if ( present(trans) .and. trans .eq. 'T' ) then
            do i = cols(1), cols(2)
                call prtVec( adns(rows(1):rows(2),i), fmt )
            end do
        else
            if ( present(fmt) ) then
                do i = rows(1), rows(2)
                    call prtVec( adns(i,cols(1):cols(2)), fmt )
                end do
            else
                write(*,*) adns
            end if
        end if
        if ( allocated(adns) ) deallocate(adns)
    end subroutine

    subroutine mso_sparse_showR(a, ia, ja, rows, cols, fmt, trans)
        use common_functions, only: prtVecR
        implicit none
        real(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        integer, intent(in) :: rows(2), cols(2)
        character(len=*), intent(in), optional :: fmt
        character(len=1), intent(in), optional :: trans

        integer :: ii, jj, mn(2), m, job(6), info, n
        real(kind=8), dimension(:,:), allocatable :: adns

        m = size(ia) - 1
        n = maxval(ja)
        allocate(adns(m,n))

        job = 1
        job(4) = 2
        call mkl_ddnscsr(job,m,n,adns,m,a,ja,ia,info)
        if ( .not.present(trans) .or. (trans.eq.'N' .or. trans.eq.'n') ) then
            do ii = rows(1), rows(2), 1
                do jj = cols(1), cols(2), 1
                    write(*,'('//fmt//'$)') adns(ii,jj)
                end do
                write(*,*)
            end do
        else
            if ( present(trans) .and. (trans.eq.'T' .or. trans.eq.'t') ) then
                do ii = cols(1), cols(2), 1
                    do jj = rows(1), rows(2), 1
                        write(*,'('//fmt//'$)') adns(jj,ii)
                    end do
                    write(*,*)
                end do
            else
                write(*,*) 'trans is not well defined in mso_sparse_showR.'
            end if
        end if
        if ( allocated(adns) ) deallocate(adns)
    end subroutine mso_sparse_showR




    subroutine mso_dgmres_nonprecondition(maxiter, maxrelres, a,ia,ja,rhs, computed_solution, itercount, res,msglvl)
        !*******************************************************************************
        !  Content:
        !  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
        !  msglvl: print the procesure data if msglvl>0
        !  res(2): res(1) is the Euclidean norm of the current residual
        !          res(2) is equal to res(1) / norm2(rhs)
        !*******************************************************************************
        use blas95, only: nrm2
        implicit none
        integer, intent(in) :: maxiter(2)
        real(kind=8), intent(in) :: maxrelres
        real(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        real(kind=8), dimension(:), intent(in) :: rhs
        real(kind=8), dimension(:), intent(inout) :: computed_solution
        integer, intent(inout) :: itercount
        real(kind=8), intent(inout) :: res(2)
        integer, optional :: msglvl

        integer :: n
        integer :: ipar(128)
        real(kind=8) :: dpar(128)
        real(kind=8), dimension(:), allocatable :: tmp
        real(kind=8), dimension(:), allocatable :: residual

        integer :: rci_request, itmp

        n = size(ia) - 1

        allocate( residual(n))

        if ( maxiter(2).gt.0 ) then
            itmp = maxiter(2)
        else
            itmp = minval((/150,n/))
        end if

        allocate( tmp( (2*itmp+1)*n+itmp*(itmp+9)/2+1 ) )

        call dfgmres_init(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999
        !---------------------------------------------------------------------------
        ! set the desired parameters:
        ! logical parameters:
        ! do residual stopping test
        ! do not request for the user defined stopping test
        ! do the check of the norm of the next generated vector automatically
        ! double precision parameters
        !---------------------------------------------------------------------------
        ipar(9) = 1
        ipar(10) = 0
        ipar(12) = 1
        if ( maxiter(1).gt.0 ) then
            ipar(5) = maxiter(1)
        else
            ipar(5) = minval((/150,n/))
        end if
        if ( maxiter(2).gt.0 ) then
            ipar(15) = maxiter(2)
        else
            ipar(15) = minval((/150,n/))
        end if
        dpar(1) = maxrelres
        !---------------------------------------------------------------------------
        ! check the correctness and consistency of the newly set parameters
        !---------------------------------------------------------------------------
        call dfgmres_check(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999
        !---------------------------------------------------------------------------
        ! print the info about the rci fgmres method
        !---------------------------------------------------------------------------
        if ( present(msglvl).and.msglvl.gt.0 ) then
            print *, ''
            print *,'some info about the current run of rci fgmres method:'
            print *, ''
            if (ipar(8) .ne. 0) then
                write(*,'(a,i1,a)') 'as ipar(8)=',ipar(8),', the automatic test for the maximal number of iterations will be performed'
                write(*,*) '(default) iteration number is ipar(5),ipar(15)', ipar(5), ipar(15)
            else
               write(*,'(a,i1,a)') 'as ipar(8)=',ipar(8),', the automatic test for the maximal number of iterations will be skipped'
            end if
            print *,'+++'
            if (ipar(9) .ne. 0) then
                write(*,'(a,i1,a)') 'as ipar(9)=',ipar(9),', the automatic residual test will be performed'
            else
                write(*,'(a,i1,a)') 'as ipar(9)=',ipar(9),', the automatic residual test will be skipped'
            end if
            print *,'+++'
            if (ipar(10) .ne. 0) then
                write(*,'(a,i1,a)') 'as ipar(10)=',ipar(10),', the user-defined stopping test will be requested via rci_request=2'
            else
                write(*,'(a,i1,a)') 'as ipar(10)=',ipar(10),', the user-defined stopping test will not be requested, thus, rci_request will not take the value 2'
            end if
            print *,'+++'
            if (ipar(11) .ne. 0) then
                write(*,'(a,i1,a)') 'as ipar(11)=',ipar(11),', the preconditioned fgmres iterations will be performed, thus, the preconditioner action will be',' requested via rci_request=3'
            else
                write(*,'(a,i1,a)') 'as ipar(11)=',ipar(11),', the preconditioned fgmres iterations will not be performed, thus, rci_request will not take the',' value 3'
            end if
            print *,'+++'
            if (ipar(12) .ne. 0) then
                write(*,'(a,i1,a)') 'as ipar(12)=',ipar(12),', the automatic test for the norm of the next generated vector is not equal to zero up to rounding and',' computational errors will be performed, thus, rci_request will not take the value 4'
            else
                write(*,'(a,i1,a)') 'as ipar(12)=',ipar(12),', the automatic test for the norm of the next generated vector is not equal to zero up to rounding and',' computational errors will be skipped, thus, the user-defined test will be',' requested via rci_request=4'
            end if
            print *,'+++'
        end if

        !---------------------------------------------------------------------------
        ! compute the solution by rci (p)fgmres solver without preconditioning
        ! reverse communication starts here
        !---------------------------------------------------------------------------
        1 call dfgmres(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
        !---------------------------------------------------------------------------
        ! if rci_request=0, then the solution was found with the required precision
        !---------------------------------------------------------------------------
        if (rci_request .eq. 0) go to 3
        !---------------------------------------------------------------------------
        ! if rci_request=1, then compute the vector a*tmp(ipar(22))
        ! and put the result in vector tmp(ipar(23))
        !---------------------------------------------------------------------------
        if (rci_request .eq. 1) then
            call mkl_dcsrgemv('n',n, a, ia, ja, tmp(ipar(22)), tmp(ipar(23)))
            go to 1
        !---------------------------------------------------------------------------
        ! if rci_request=anything else, then dfgmres subroutine failed
        ! to compute the solution vector: computed_solution(n)
        !---------------------------------------------------------------------------
        else
            go to 999
        end if
        !---------------------------------------------------------------------------
        ! reverse communication ends here
        ! get the current iteration number and the fgmres solution (do not forget to
        ! call dfgmres_get routine as computed_solution is still containing
        ! the initial guess!)
        !---------------------------------------------------------------------------
        3 call dfgmres_get(n, computed_solution, rhs, rci_request, ipar, dpar, tmp, itercount)
        call mkl_dcsrgemv('n', n, a, ia, ja, computed_solution, residual)
        call daxpy(n, -1.0d0, rhs, 1, residual, 1)
        res(1) = nrm2(residual)
        res(2) = res(1)/nrm2(rhs)
        !---------------------------------------------------------------------------
        ! print the number of iterations: itercount
        !---------------------------------------------------------------------------
        if ( present(msglvl).and.msglvl.gt.0 ) then
            print *, ''
            print *,' the system has been solved'
            print *, ''
            print *,' number of iterations: ', itercount
            print *,'the absolute tolerance is ', res(1)
            print *,'the relative tolerance is ', res(2)
        end if
        if ( allocated(residual) ) deallocate(residual)
        if ( allocated(tmp) ) deallocate(tmp)
        !---------------------------------------------------------------------------
        ! release internal intel(r) mkl memory that might be used for computations
        ! note: it is important to call the routine below to avoid memory leaks
        ! unless you disable intel(r) mkl memory manager
        !---------------------------------------------------------------------------
        call mkl_free_buffers
        return

        !---------------------------------------------------------------------------
        ! release internal intel(r) mkl memory that might be used for computations
        ! note: it is important to call the routine below to avoid memory leaks
        ! unless you disable intel(r) mkl memory manager
        !---------------------------------------------------------------------------
        999 if ( rci_request.eq.-1 ) then
            call dfgmres_get(n, computed_solution, rhs, rci_request, ipar, dpar, tmp, itercount)
            call mkl_dcsrgemv('n', n, a, ia, ja, computed_solution, residual)
            call daxpy(n, -1.0d0, rhs, 1, residual, 1)
            res(1) = nrm2(residual)
            res(2) = res(1)/nrm2(rhs)

            write( *,'(a,i5)') 'this example failed as the solver has returned the error code', rci_request
            print *, ''
            print *,' the system has not been solved because the maximum number of iterations was reached, but the relative stopping criterion was not met'
            print *, ''
            print *,'current number of iterations: ', itercount
            print *,'the absolute tolerance is ', res(1)
            print *,'the relative tolerance is ', res(2)
            if ( allocated(residual) ) deallocate(residual)
            if ( allocated(tmp) ) deallocate(tmp)
            call mkl_free_buffers
            return
        else
            write( *,'(a,i5)') 'this example failed as the solver has returned the error code', rci_request
            write(*,*) 'some errors accurred in mso_dgmres_nonprecondition!'
            if ( allocated(residual) ) deallocate(residual)
            if ( allocated(tmp) ) deallocate(tmp)
            call mkl_free_buffers
            stop 1
        end if
    end subroutine mso_dgmres_nonprecondition

    subroutine mso_dgmres_precondition(maxiter, relres, a,ia,ja,rhs, computed_solution, itercount)
        !*******************************************************************************
        !  Content:
        !  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
        !*******************************************************************************
        implicit none
        integer, intent(in) :: maxiter
        real(kind=8) :: relres
        real(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        real(kind=8), dimension(:), intent(in) :: rhs
        real(kind=8), dimension(:), intent(inout) :: computed_solution
        integer, intent(in) :: itercount


        integer :: n, m
        integer :: ipar(128)
        real(kind=8) :: dpar(128)
        real(kind=8), dimension(:), allocatable :: tmp, b, residual

        integer :: rci_request, i
        real(kind=8) :: dvar
        real(kind=8), external :: dnrm2
        real(kind=8), dimension(:), allocatable :: bilu0, trvec
        integer :: ierr, err

        n = size(ia) - 1    ! the size of the problem

        allocate( b(n))
        allocate( residual(n))

        !---------------------------------------------------------------------------
        ! Initialize variables and the right hand side through matrix-vector product
        !---------------------------------------------------------------------------

        call dcopy(n, rhs, 1, b, 1)

        computed_solution = 1.d0

        ipar(15) = min(150,n,maxiter)
        allocate( tmp(n*(2*ipar(15)+1)+(ipar(15)*(ipar(15)+9))/2+1))

        call dfgmres_init(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        ipar(8) = 1 ! if the value is not equal to 0, the dfgmres routine performs the stopping test for the maximum number of iterations: ipar(4) ≤ ipar(5)
        ipar(15) = min(150,n,maxiter)
        ipar(5) = ipar(15)
        ipar(11) = 0
        dpar(1) = relres

        call dfgmres_check(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        ipar(31) = 1
        dpar(31) = 1.D-20
        dpar(32) = 1.D-16
        allocate(bilu0(size(a)), stat=err)
        if (err /= 0) print *, "bilu0: Allocation request denied"
        call dcsrilu0(n, a, ia, ja, bilu0, ipar, dpar, ierr)
        allocate(trvec(n))

        if (ierr .ne. 0) then
            write(*,'(a,a,i1)') ' error after calculation of the',' preconditioner dcsrilu0',ierr
              go to 999
        end if

        1 call dfgmres(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)

        if (rci_request .eq.0) go to 3

        if (rci_request .eq. 1) then
            call mkl_dcsrgemv('n',n, a, ia, ja, tmp(ipar(22)),tmp(ipar(23)))
            go to 1
        end if

        if (rci_request .eq. 2) then
            ! request to the dfgmres_get routine to put the solution into b(n) via ipar(13)
            ipar(13) = 1
            ! get the current fgmres solution in the vector b(n)
            call dfgmres_get(n, computed_solution, b, rci_request, ipar,dpar, tmp, itercount)
            ! compute the current true residual via mkl (sparse) blas routines
            call mkl_dcsrgemv('n', n, a, ia, ja, b, residual)
            call daxpy(n, -1.0d0, rhs, 1, residual, 1)
            dvar = dnrm2(n, residual, 1)
            if (dvar .lt. relres) then
                go to 3
            else
                go to 1
            end if
        end if

        if (rci_request .eq. 3) then
            call mkl_dcsrtrsv('l','n','u',n,bilu0,ia,ja,tmp(ipar(22)),trvec)
            call mkl_dcsrtrsv('u','n','n',n,bilu0,ia,ja,trvec,tmp(ipar(23)))
            go to 1
        end if

        if (rci_request .eq. 4) then
            if (dpar(7) .lt. 1.0d-12) then
                go to 3
            else
                go to 1
            end if
            !---------------------------------------------------------------------------
            ! if rci_request=anything else, then dfgmres subroutine failed
            ! to compute the solution vector: computed_solution(n)
            !---------------------------------------------------------------------------
        else
            go to 999
        end if

        3 IPAR(13) = 0
        call dfgmres_get(n, computed_solution, rhs, rci_request, ipar, dpar, tmp, itercount)


        if ( allocated(tmp) ) deallocate(tmp)
        if ( allocated(bilu0) ) deallocate(bilu0)
        if ( allocated(trvec) ) deallocate(trvec)
        if ( allocated(b) ) deallocate(b)
        if ( allocated(residual) ) deallocate(residual)
        return

        999 write( *,'(A,I5)') 'this example failed as the solver has returned the error code', rci_request
        call mkl_free_buffers
        stop 1
    end subroutine mso_dgmres_precondition

    subroutine mso_dpardiso(mtype,msglvl,a,ia,ja,b,x)
        ! mtype: 2 for real(8) and symmetric positive definite
        !        -2 for real(8) and symmetric indefinite
        !        11 for real(8) and nonsymmetric
        ! msglvl: 0 for no information output
        !         1 for information output
        USE mkl_pardiso
        use common_functions, only: cf_allocate, cf_deallocate
        IMPLICIT NONE
        integer, parameter :: dp = kind(1.d0)
        integer, intent(in) :: mtype, msglvl
        integer, dimension(:), intent(in) :: ia, ja
        real(kind=dp), dimension(:), intent(in) :: a
        real(kind=dp), dimension(:), intent(inout) :: b
        real(kind=dp), dimension(:), intent(inout) :: x

        !.. internal solver memory pointer
        type(mkl_pardiso_handle), allocatable :: pt(:)
        !.. all other variables
        integer :: maxfct, mnum, phase, n, nrhs, error, nnz
        integer :: error1
        integer, allocatable :: iparm(:)
        integer :: i, idum(1)
        real(kind=dp) :: ddum(1)
        !.. Fill all arrays containing matrix data.
        n = size(ia) - 1
        nnz = size(a)
        nrhs = 1
        maxfct = 1
        mnum = 1

        allocate(pt(64))
        call cf_allocate(iparm,64)
        !..
        !.. set up pardiso control parameter
        !..
        !         allocate(iparm(64))

        do i = 1, 64
           iparm(i) = 0
        end do

        iparm(1) = 1 ! no solver default
        iparm(2) = 2 ! fill-in reordering from METIS
        iparm(4) = 0 ! no iterative-direct algorithm
        iparm(5) = 0 ! no user fill-in reducing permutation
        iparm(6) = 0 ! =0 solution on the first n components of x
        iparm(8) = 2 ! numbers of iterative refinement steps
        iparm(10) = 13 ! perturb the pivot elements with 1E-13
        iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
        iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
        iparm(14) = 0 ! Output: number of perturbed pivots
        iparm(18) = -1 ! Output: number of nonzeros in the factor LU
        iparm(19) = -1 ! output: mflops for lu factorization
        iparm(20) = 0 ! output: numbers of cg iterations

        error  = 0 ! initialize error flag
        !       msglvl = 1 ! print statistical information
        !       mtype  = -2 ! symmetric, indefinite

        !.. initialize the internal solver memory pointer. this is only
        ! necessary for the first call of the pardiso solver.

        !         allocate (pt(64))
        do i = 1, 64
          pt(i)%dummy =  0
        end do

        !.. back substitution and iterative refinement
        iparm(8) = 0 ! max numbers of iterative refinement steps
        phase = 13 ! 33 ! only solving
        write(*,*) 'pardiso start'
        call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
              idum, nrhs, iparm, msglvl, b, x, error)
        ! write(*,*) 'solve completed ... '
        if (error /= 0) then
            write(*,*) 'the following error was detected: ', error
            goto 1000
        end if

        1000 continue
        write(*,*) '10000000'
        !.. termination and release of memory
        phase = -1 ! release internal memory
        call pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
              idum, nrhs, iparm, msglvl, ddum, ddum, error1)

        !         if (allocated(iparm))   deallocate(iparm)
        deallocate(pt)
        call cf_deallocate(iparm)
        if (error1 /= 0) then
            write(*,*) 'the following error on release stage was detected: ', error1
            stop 1
        endif

        if (error /= 0) stop 1
    end subroutine mso_dpardiso

    subroutine mso_zpardiso(mtype,msglvl,a,ia,ja,b,x)
        USE mkl_pardiso
        IMPLICIT NONE
        integer, intent(in) :: mtype, msglvl
        integer, dimension(:), intent(in) :: ia, ja
        complex(kind=8), dimension(:), intent(in) :: a
        complex(kind=8), dimension(:), intent(inout) :: b
        complex(kind=8), dimension(:), intent(inout) :: x

        !.. internal solver memory pointer
        type(mkl_pardiso_handle), allocatable  :: pt(:)
        !.. all other variables
        integer :: maxfct, mnum, phase, n, nrhs, error
        integer :: error1
        integer, allocatable :: iparm( : )
        integer :: i, idum(1)
        real(kind=8) :: ddum(1)
        !.. Fill all arrays containing matrix data.
        n = size(ia) - 1
        nrhs = 1
        maxfct = 1
        mnum = 1
        !..
        !.. set up pardiso control parameter
        !..
        allocate(iparm(64))

        do i = 1, 64
           iparm(i) = 0
        end do

        iparm(1) = 1 ! no solver default
        iparm(2) = 2 ! fill-in reordering from METIS
        iparm(4) = 0 ! no iterative-direct algorithm
        iparm(5) = 0 ! no user fill-in reducing permutation
        iparm(6) = 0 ! =0 solution on the first n components of x
        iparm(8) = 2 ! numbers of iterative refinement steps
        iparm(10) = 13 ! perturb the pivot elements with 1E-13
        iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
        iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
        iparm(14) = 0 ! Output: number of perturbed pivots
        iparm(18) = -1 ! Output: number of nonzeros in the factor LU
        iparm(19) = -1 ! output: mflops for lu factorization
        iparm(20) = 0 ! output: numbers of cg iterations

        error  = 0 ! initialize error flag
        !       msglvl = 1 ! print statistical information
        !       mtype  = -2 ! symmetric, indefinite

        !.. initialize the internal solver memory pointer. this is only
        ! necessary for the first call of the pardiso solver.

        allocate (pt(64))
        do i = 1, 64
          pt(i)%dummy =  0
        end do

        !.. back substitution and iterative refinement
        iparm(8) = 2 ! max numbers of iterative refinement steps
        phase = 13 ! 33 ! only solving

        call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, b, x, error)
        ! write(*,*) 'solve completed ... '
        if (error /= 0) then
            write(*,*) 'the following error was detected: ', error
            goto 1000
        end if

        1000 continue
        !.. termination and release of memory
        phase = -1 ! release internal memory
        call pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
              idum, nrhs, iparm, msglvl, ddum, ddum, error1)

        if (allocated(iparm))   deallocate(iparm)
        if (allocated(pt))   deallocate(pt)

        if (error1 /= 0) then
            write(*,*) 'the following error on release stage was detected: ', error1
            stop 1
        endif

        if (error /= 0) stop 1
    end subroutine mso_zpardiso

end module MySparseOperator